#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018/1/12 13:30
# @Author  : Deyu
# @Site    : 
# @File    : explory_data.py
# @Software: PyCharm Community Edition
# @Function explorery data analysis

from __future__ import print_function, division
import math
import numpy as np
import numpy.ma as ma
import scipy

from pandas import  DataFrame as DF
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy.stats import pearsonr

from geotiff_tools import *
from config import *

# add folder to the Python path at runtime
#从其他文件夹导入package
import sys
#sys.path.append('''P:\\Aoi_paper\\codes\\ThinkStats2\\code''')
# sys.path.append('''..\\ThinkStats2\\code''')
#print(sys.path)
# import thinkstats2 as thst
# from thinkstats2 import thinkplot as thplt



def stackNumricFeatures1(numric, imggt):
    """
    合并类jian相关性检测后保留的9个连续变量
    :param basic:
    :param hydro:
    :param meropy:
    :param geos:
    :return:
    """
    pass
    print(numric.shape)
    tmp = np.stack((numric[0], numric[1], numric[2], numric[3], numric[4], numric[5], numric[6], numric[7], numric[11]))
    print(tmp.shape)
    for i in range(len(tmp)):
        array2rasterUTM("{}/{}-numric-features.tif".format(config.outdir, i+1), imggt, tmp[i])



def interClassPearsonRsquareTest(numric_feature):
    """
    r2 correlection interblock
    :param basic:
    :param hydro:
    :param meropy:
    :param geos:
    :return:
    """
    for i in range(13):
        for j in range(13):
            print(i, j, pearsonr(numric_feature[i].ravel(), numric_feature[j].ravel())[0])
    pass





def stackNumricFeatures(basic, hydro, meropy, geos, imggt):
    """
    合并类内相关性检测后保留的13个连续变量
    :param basic:
    :param hydro:
    :param meropy:
    :param geos:
    :return:
    """
    pass
    print(meropy.shape, geos.shape, hydro.shape, basic.shape)
    tmp = np.stack((basic[0], basic[1], basic[2], basic[4], hydro[0], hydro[3], meropy[0], meropy[1], meropy[2], meropy[3], geos[0], geos[2], geos[4]))
    print(tmp[0].shape)
    for i in range(len(tmp)):
        array2rasterUTM("{}/{}-numric-features.tif".format(config.outdir, i+1), imggt, tmp[i])

def inClassPearsonRsquareTest(basic, hydro, meropy, geos):
    """
    r2 correlection in class
    :param basic:
    :param hydro:
    :param meropy:
    :param geos:
    :return:
    """
    for i in range(5):
        for j in range(5):
            print(i, j, pearsonr(basic[i].ravel(), basic[j].ravel())[0])
    pass

    for i in range(4):
        for j in range(4):
            print(i, j, pearsonr(hydro[i].ravel(), hydro[j].ravel())[0])
    pass

    for i in range(4):
        for j in range(4):
            print(i, j, pearsonr(meropy[i].ravel(), meropy[j].ravel()))
    pass

    for i in range(5):
        for j in range(5):
            print(i, j, pearsonr(geos[i].ravel(), geos[j].ravel()))
    pass




def PearsonChisquareTest(mory, geos):
    """
    chi-square test for nominal feature correlation
    :param arr:
    :return:
    """
    TPI = mory[4].ravel()
    IP = mory[5].ravel()
    FOS = geos[1].ravel()

    arr1 = np.stack((IP, FOS), axis=1)
    print(arr1.shape)
    df = DF(arr1, columns=('IP', 'FOS'), dtype=np.int64)
    df_samp = df.sample(n=300)
    table = sm.stats.Table.from_data(df_samp[["IP", "FOS"]])
    print(table.table_orig)
    print('statistic', table.test_nominal_association().statistic)
    print('degree of freedom', table.test_nominal_association().df)
    print('pvalue', table.test_nominal_association().pvalue)
    pass
    #print('chi-square-0.95-28:',scipy.stats.chi2.isf(0, 28))
    #print('chi-square of p value with 1 and degree of freedom 36:', scipy.stats.chi2.isf(0, 36))
def plot_PMF_pyplot():
    """
    :param:
    坡度和坡向的概率质量函数
    :return:
    """
    arr = read_tiff("{}/slope_aspect_curvature_vallydepth.tif".format(terrainDir))[0]
    arr1 = read_tiff("{}/slope_aspect_curvature_vallydepth.tif".format(terrainDir))[1]
    x = arr.ravel()
    x = x[x > 1e-4]
    x = x * 180 / math.pi
    print(np.max(x), np.min(x))

    y = arr1.ravel()
    y = y[y > -1]
    #print(len(np.unique(y))) #1940
    y = y * 180 / math.pi
    print(np.max(y), np.min(y))

    # the histogram of the data
    n, bins, patches = plt.hist(x, bins=200, normed=1,
                                linewidth=2, color='tab:red', alpha=0.8, label='slope',
                                stacked=1, cumulative=0, histtype='bar')

    n, bins, patches = plt.hist(y, bins=200, normed=1,
                                linewidth=2, color='m', alpha=0.7, label='aspect',
                                stacked=1, cumulative=0, histtype='bar')
    print(n)  #print(bins)  #print(patches)
    plt.xlabel('Surface(Degree)', fontsize=10.5)
    plt.ylabel('Probability', fontsize=10.5)
    plt.xticks(fontsize=10.5)
    plt.yticks(np.arange(0, 0.05, 0.01), fontsize=10.5)
    plt.xticks(np.arange(0, 360, 90), fontsize=10.5)
    plt.axis([0, np.max(y), 0, 0.05])
    plt.legend(prop={'size': 10.5}, loc='upper right')
    plt.grid(True)

    plt.show()
    pass
def plot_PMF_thinkstats2():
    arr = img2array("{}\\table_4_1_basic_terrain.tif".format(config.terrainDir))
    print(arr[0][0].shape)
    slope = arr[0].ravel()
    #slope =arr[0][1000]
    #slope[slope<0] = 0
    slope = ma.masked_where(slope < 0, slope)
    slope = np.asarray(slope)
    print(np.max(slope), np.min(slope))
    pmf = thst.Pmf(slope)
    print(pmf.Total())
    pmf.Normalize()
    thplt.PrePlot(1)
    thplt.Hist(pmf)
    thplt.Config(xlabel='slope', ylabel='Prob Mass Func')
    thplt.show()

    pass


if __name__ == '__main__':
    # imggt = read_tif_metadata("{}/stack-numric-features.tif".format(config.terrainDir))
    # basic_terrain = read_tiff("{}/table_4_1_basic_terrain.tif".format(config.terrainDir))
    # hydro_terrain = read_tiff("{}/table_4_2_hydro_terrain.tif".format(config.terrainDir))
    # meropy_terrain = read_tiff("{}/table_4_3_morpy_terrain.tif".format(config.terrainDir))
    # geos_terrain = read_tiff("{}/table_4_5_hygeos_terrain.tif".format(config.terrainDir))
    # numric_features = read_tiff("{}/stack-numric-features.tif".format(config.terrainDir))
    # print(np.unique(meropy_terrain[4]))
    # print(np.unique(meropy_terrain[5]))
    # print(np.unique(geos_terrain[1]))
    #PearsonChisquareTest(meropy_terrain, geos_terrain)
    #plot_PMF_pyplot(basic_terrain)
    #plot_PMF_thinkstats2()
    #inClassPearsonRsquareTest(basic_terrain, hydro_terrain, meropy_terrain, geos_terrain)
    #stackNumricFeatures(basic_terrain, hydro_terrain, meropy_terrain, geos_terrain, imggt)
    #interClassPearsonRsquareTest(numric_features)
    # stackNumricFeatures1(numric_features, imggt)
    plot_PMF_pyplot()
    pass

